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Abstract 


This  report  describes  the  temperature  distribution  within  the  Small-scale  Cookoff  Bomb 
(SCB)  test  apparatus  from  both  an  experimental  and  theoretical  perspective.  The  two- 
dimensional  finite  difference  computer  code  HEAT  is  used  to  predict  temperature-time 
profiles  at  specified  positions  within  an  inert  filled  SCB.  These  predictions  are  then 
compared  with  temperatures  obtained  from  thermocouple  measurements  made  on  the 
test  apparatus.  The  general  agreement  is  good  and  justifies  the  continued  application 
of  HEAT  to  the  SCB.  Some  suggestions  for  further  modifications  to  HEAT  are  made 
which  would  increase  the  accuracy  of  the  calculation. 
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Heat  Flow  Calculations  for 
the  Small-Scale  Cookoff  Bomb  Test 


1.  Introduction 

A  Super  Small-scale  Cookoff  Bomb  (SSCB)  test  facility  has  recently  been 
established  at  MRL  to  provide  data  on  the  cookoff  behaviour  of  energetic 
materials  [1].  This  test  is  suitable  for  explosive  samples  having  a  mass  of  about 
20  g.  The  behaviour  of  explosive  samples  of  considerably  larger  mass  (up  to 
700  g)  is  evaluated  using  the  Small-scale  Cookoff  Bomb  (SCB)  test  [2].  The 
development  of  an  SCB  test  at  MRL  is  well  advanced. 

As  part  of  the  program  to  assess  the  cookoff  behaviour  of  both  existing  and 
new  explosive  formulations  at  MRL  we  have  made  some  preliminary  heat  flow 
calculations  for  the  MRL  SCB  test  apparatus.  The  calculations  were  performed 
with  the  HEAT  computer  code  [31  to  ascertain  the  suitability  of  the  code  for 
modelling  the  thermal  environment  in  the  SCB  test. 

HEAT  is  a  reactive  multi-material  finite  difference  Fortran  computer  code  for 
the  solution  cf  the  diffusion  equation  in  two-dimensional  axisymmetric 
cylindrical  geometry.  It  allows  a  maximum  of  three  separate  nested  materials, 
and  in  a  typical  application  the  innermost  material  will  be  an  explosive.  HEAT 
therefore  allows  the  innermost  material  to  undergo  a  change  in  phase,  and  also 
to  decompose  via  zero  order  Arrhenius  kinetics. 

The  design  of  the  SCB  is  shown  in  Figure  1.  It  is  constructed  from  a  single 
steel  cylinder  with  a  welded  steel  base.  The  upper  section  of  the  cylinder  is 
threaded  to  allow  it  to  be  sealed  by  a  screw-on  steel  cap.  The  SCB  is  then 
placed  between  two  127  mm  square,  12.7  mm  thick  steel  plates  which  are 
secured  by  four  bolts.  A  small  air  gap  exists  between  the  upper  surface  of  the 
explosive  and  the  steel  cap.  Numerical  modelling  of  the  temperature 
distribution  within  the  SCB  therefore  requires  a  computer  code  with  the 
capability  of  simulating  heat  flow  through  the  steel  cylinder,  the  explosive,  and 
the  air  gap.  HEAT  can  model  three  nested  materials  and  can  be  adapted  to  the 
geometry  shown  in  Figure  1,  and  was  therefore  used  for  our  numerical 
modelling  of  heat  transfer  in  cookoff  bomb  tests  by  simulating  the  heat  flow  in 
the  SCB.  Numerical  simulation  of  heat  flow  in  the  SSCB  will  be  discussed  in 
a  future  report. 
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Figure  1:  Small-scale  Cookoff  Bomb  (SCB). 


HEAT  was  originally  developed  to  study  the  transient  heating  of  a  shell 
within  a  hot  gun  barrel.  It  employs  the  cylindrical  geometry  of  a  shell,  with 
the  boundary  conditions  imposing  a  time-independent  temperature  profile  over 
the  top  and  bottom  surfaces  of  the  outermost  material.  A  time  dependent 
programmed  temperature  is  allowed  along  the  curved  cylindrical  surface. 
Preliminary  experimentation  with  the  SCB  using  an  inert  explosive  replacement 
and  typical  heating  rates  showed  that  the  SCB  behaviour  was  not  compatible 
with  the  assumption  of  a  time  independent  temperature  profile  over  the  end 
surfaces.  Heat  loss  over  these  surfaces  from  both  radiation  and  convection  was 
also  noted.  Changes  to  the  code  were  therefore  made  to  alter  the  boundary 
conditions  over  these  surfaces  to  allow  a  time  dependent  variable  temperature 
profile  and  heat  losses  due  to  both  radiation  and  convection.  Because  of  the 
method  used  to  solve  the  finite  difference  equations  and  the  nonlinear  nature  of 
the  radiation  boundary  condition,  this  entailed  making  quite  detailed  changes  to 
the  code.  These  changes  are  described  in  Section  2,  along  with  the  procedure 
we  used  to  estimate  realistic  values  for  both  the  emissivity  and  heat  transfer 
coefficient  for  the  steel  surfaces  of  the  SCB. 

To  provide  an  experimental  data  base  for  comparison  with  the  code 
calculations  a  series  of  SCB  test  runs  was  made  using  sand  as  an  inert  explosive 
substitute.  Two  Type  K  thermocouples  were  used  in  each  test  to  record  the 
temperature  at  various  locations  within  the  SCB;  typically  these  were  at  the 


sand /steel  interface,  and  either  the  centre  of  the  sand,  or  the  surface  of  the  steel 
cap.  These  locations  are  indicated  on  the  schematic  of  the  SCB  in  Figure  1. 

The  thermocouples  were  connected  to  dedicated  thermocouple  amplifiers  and 
the  output  from  these  amplifiers  was  then  recorded  using  a  strip  chart  recorder. 

In  order  to  obtain  accurate  estimates  of  the  surface  heat  transfer  coefficient 
and  emissivity  of  the  SCB  surface  a  further  series  of  test  runs  was  made  on  a 
solid  steel  cylinder  having  approximately  the  same  dimensions  as  the  SCB.  Up 
to  nine  thermocouples  were  used  in  each  of  the  runs  to  measure  temperature  at 
various  positions  and  depths  into  the  cylinder,  and  the  temperatures  were 
recorded  using  a  DT  100  Datataker. 

In  Section  4  a  comparison  is  made  between  the  HEAT  code  simulations  and 
selected  experimental  records  from  the  above  two  series  of  runs.  In  general  the 
agreement  is  good,  and  justifies  the  continued  application  of  HEAT  to  the 
numerical  simulation  of  the  SCB.  The  agreement  is  not  perfect,  however,  and  in 
most  cases  there  are  differences  of  the  order  of  10°C  in  some  sections  of  the 
temperature-time  records.  A  possible  reason  for  this  discrepancy  is  the  neglect 
of  the  temperature  dependence  of  the  thermal  diffusivity  of  each  of  the 
materials.  This  is  also  discussed  further  in  Section  4.  In  Section  5  we 
summarize  our  results,  and  discuss  directions  for  future  work. 


2.  HEAT  Code  Modifications 

The  HEAT  code  has  previously  been  described  in  detail  [3],  It  is  a  finite 
difference  code  which  solves  the  partial  differential  equation  describing  time 
dependent  thermal  conduction  with  an  internal  heat  source,  i.e. 

?I--L.V2T  +  S  (1) 

7T  JF 


where  T  is  the  temperature,  p  the  density,  X  the  thermal  conductivity,  and  Cp 
the  specific  heat  at  constant  pressure.  S  is  the  source  term  representing 
internal  heat  generation  due  to  thermal  decomposition  of  the  explosive.  In  the 
calculations  described  here  we  consider  only  inert  materials  and  so  have  set  S 
equal  to  zero.  The  combination  X/(,pCf)  is  called  the  thermal  diffusivity  of  the 
material.  The  value  of  this  parameter  determines  the  time  dependence  of  the 
temperature  distribution  within  the  material. 

Equation  (1)  is  solved  in  two-dimensional  axisymmetric  cylindrical  geometry 
with  r  and  z  representing  displacements  in  the  radial  and  axial  directions 
respectively.  A  second  order  Forward-Time  Centred-Space  (FTCS)  scheme  is 
used  to  discretize  the  equations  [4],  and  then  the  resulting  set  of  algebraic 
equations  is  solved  using  the  Alternating  Direction  Implicit  (ADI)  method  [51. 
This  is  a  combination  of  implicit  and  explicit  methods  which  solves  the 
equations  in  two  passes  during  one  time  step.  In  the  first  pass  the  partial 
derivatives  in  one  direction  (say  z)  are  replaced  by  a  finite  difference  expression 
in  terms  of  known  values  of  T  (explicit),  while  the  remaining  partial  derivatives 
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(r  direction)  are  replaced  by  finite  difference  expressions  in  terms  of  unknown 
values  (implicit).  In  the  second  pass  the  z  direction  is  made  implicit  and  the  r 
direction  explicit.  The  ADI  method  has  several  advantages  over  alternative 
finite-difference  schemes:  (a)  the  resulting  sets  of  equations  have  a  simple 
tridiagonal  matrix  form  and  so  can  easily  be  solved  and  (b)  the  equations  are 
reported  [4-6]  to  have  stable  solutions  for  any  combination  of 
time  step  size  and  spatial  cell  size. 

In  our  early  experimentation  with  HEAT  we  found  that  it  was  possible  to 
choose  combinations  of  At,  Ar  and  Az  which  gave  unstable  solutions  for  the 
temperature  profiles  and  led  to  negative  values  for  the  absolute  temperature. 

In  all  cases  we  were  able  to  trace  the  origin  of  these  unstable  solutions  to  the 
presence  of  boundaries  between  different  materials.  While  the  ADI  method  is 
known  to  be  unconditionally  stable  within  a  given  material  it  is  perhaps  not 
surprising  that  this  stability  breaks  down  at  an  interface  because  the  finite 
difference  schemes  used  in  HEAT  are  second  order  within  the  material,  but 
only  first  order  at  the  interface.  The  original  report  [3]  noted  this  problem  and 
gave  guidelines  for  the  avoidance  of  the  unstable  solutions.  In  particular,  it 
was  recommended  that  at  least  10  computational  cells  span  a  given  material  in 
any  one  direction.  This  has  forced  us  to  use  a  very  fine  finite  difference  grid  to 
model  the  SCB,  but  by  complying  with  this  guideline,  and  using  realistic  values 
for  Af,  we  have  avoided  the  occurrence  of  unstable  solutions. 

The  boundary  conditions  in  the  original  version  of  HEAT  consisted  of  a  fixed 
temperature  along  the  top  and  bottom  of  the  finite  difference  grid  and  a 
prescribed  time-dependent  temperature  over  the  outer  curved  surface.  While 
this  was  adequate  for  the  original  problems  to  which  HEAT  was  applied  it  was 
considered  unsuitable  for  modelling  either  the  SCB  or  SSCB  as  it  made  no 
allowance  for  thermal  losses  due  to  radiation,  conduction  or  convection  from 
the  top  and  bottom  surfaces  of  the  apparatus.  Consequently  we  have  modified 
HEAT  so  that  the  temperatures  along  the  top  and  bottom  of  the  finite  difference 
grid  are  determined  by  a  boundary  condition  of  the  form 
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where  n  is  the  outward  normal  to  the  surface.  The  surface  heat  transfer 
coefficient  is  ft,  e  is  the  emissivity  of  the  material,  a  is  the  Stefan-Boltzmann 
constant,  and  TB  and  TA  represent  the  temperature  on  the  boundary  and  the 
ambient  temperature  of  the  surrounding  medium  (air)  respectively.  The  first 
term  on  the  right  of  equation  (2)  represents  convective  heat  loss,  while  the 
second  term  represents  heat  loss  due  to  thermal  radiation. 

To  implement  equation  (2)  we  used  the  concept  of  a  fictitious  node  point  F 
outside  the  surface  of  the  SCB.  When  the  finite  difference  approximation  to 
equation  (1)  is  applied  to  a  boundary  node  B  this  automatically  introduces  the 
temperature  Tr  at  the  fictitious  point.  By  applying  a  central  difference 
approximation  to  dT/dn  at  the  boundary  node  B,  a  second  equation  is  obtained 
which  also  contains  the  temperature  at  the  fictitious  node  point.  Tf  is  then 
eliminated  between  these  two  equations  and  a  single  equation  for  the  unknown 
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T g  is  then  obtained. 

The  solution  of  the  equation  for  Tg  is  complicated  by  the  nonlinearity  of  the 
second  term  on  the  right  of  equation  (2).  We  resolved  this  problem  by 
linearizing  the  term  using  the  following  approximation, 

V  “  Vo  +  4  Vo  <Ts  ' 

=  4  Tg30TB  -  3  Vo  (3) 

where  T8f)  is  the  temperature  at  the  boundary  node  at  the  previous  time  step. 
This  procedure  is  used  and  recommended  by  Croft  and  Lilley  [61. 

Although  the  replacement  of  the  fixed  temperature  boundary  condition  by  a 
condition  on  the  gradient  of  the  temperature  at  the  boundary  as  given  by 
equation  (2)  would  appear  to  be  fairly  straightforward,  it  is  worth  emphasizing 
that  this  resulted  in  substantial  changes  to  parts  of  the  HEAT  code.  In 
particular,  all  Fortran  DO  LOOPS  along  the  axial  direction  were  altered  to  allow 
for  additional  points  along  the  top  and  bottom  of  the  grid.  In  addition,  the 
arrays  involved  in  the  tridiagonal  matrix  calculations  had  to  be  expanded  to 
allow  for  these  additional  points. 

Before  using  the  new  version  of  HEAT,  realistic  values  for  the  coefficients  h 
and  e  had  to  be  determined.  We  did  this  by  taking  a  solid  steel  cylinder  of 
axial  length  125  mm  and  diameter  70  mm,  attaching  the  SCB  heating  bands  to 
it,  and  rapidly  heating  it  to  a  temperature  of  236°C.  We  then  monitored  the 
temperature  of  the  block  using  several  thermocouples  at  various  depths  into  the 
block  and  recorded  the  temperature-time  history  until  the  block  had  almost 
reached  ambient  temperature  again.  Because  of  the  high  thermal  conductivity 
of  the  steel  and  the  comparatively  long  time  scale  of  the  cooling  curve  all  points 
within  the  block  have  the  same  temperature  at  any  given  time  (the 
thermocouple  measurements  confirmed  this).  This  temperature  is  given  by  the 
solution  of  the  following  ordinary  differential  equation. 


mCp(dT/dt)  -  -hA(T-TA) -  eo  ACT4  -  T4 ) 


(4) 


where  A  is  the  total  surface  area  of  the  cylindrical  block. 

Using  the  values  p  =  7680  kg/m3,  Cp  =  485  J/kg.K,  and  TA  =  20°C  we  solved 
equation  (4)  using  a  simple  fourth  order  Runge-Kutta  integration  routine  for 
various  realistic  estimates  of  the  coefficients  h  and  e.  Figure  2  shows  the 
solution  obtained  using  h  =  7  W/m2  K  and  e  =  0.8  and  also  the  experimental 
points.  The  agreement  is  excellent,  and  we  have  used  these  values  for  h  and  e 
in  all  the  calculations  described  in  this  report.  The  value  for  h  obtained  here  is 
also  in  excellent  agreement  with  the  value  quoted  in  [8]  for  heat  transfer  from 
the  surface  of  a  cylindrical  steel  barrel  over  a  similar  temperature  range.  It 
should  be  noted  that,  since  the  value  of  e  depends  significantly  on  the 
oxidation  state  of  the  metal  surface,  the  value  of  e  obtained  for  the  solid  steel 
block  will  not  necessarily  be  valid  for  the  surface  of  the  steel  cap  used  on  the 
SCB,  or  the  surfaces  of  the  confining  steel  end  plates.  This  is  not  crucial 
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however  as  we  noted  very  little  change  to  the  curve  shown  in  Figure  2  when 
we  varied  e  over  the  range  0.7  to  0.9.  This  simply  reflects  the  fact  that 
convection,  rather  than  radiation,  is  the  more  important  heat  loss  mechanism  in 
the  temperature  range  shown  in  Figure  2.  The  calculations  described  in  the 
next  section  will  also  show  that  heat  loss  through  the  ends  of  the  SCB  has  only 
a  minor  effect  on  the  temperature  distribution  within  the  SCB,  and  so  very 
accurate  determinations  of  £  and  h  are  not  required. 


Figure  2:  Cooling  curve  for  solid  steel  cylinder.  Solid  line  -  from  solution  of 
ordinary  differential  equation  describing  heat  loss  due  to  convection  and  radiation. 
Filled  circles  —  experimental  pain* 


3.  Experimental  Measurements 

Before  discussing  any  of  the  HEAT  code  simulations  we  first  describe  some  of 
the  experimental  results  for  the  SCB  using  sand  as  an  explosive  substitute.  The 
upper  curve  in  Figure  3a  shows  the  temperature  of  the  sand/steel  interface  at  a 
point  located  at  about  half  the  height  of  the  cylinder,  while  the  lower  curve  is 
the  temperature  on  the  top  surface  of  the  SCB  approximately  at  the  centre.  The 
Variac  potential  was  set  to  120  volts  and  the  upper  curve  shows  that  power  was 
supplied  for  about  550  seconds  and  then  switched  off.  This  heated  the 
sand/steel  interface  to  approximately  350°C,  and  then  very  rapid  cooling  of  the 
steel  cylinder  occurred.  The  temperature  on  the  upper  surface  of  the  SCB  rose 
very  slowly  however,  reaching  a  temperature  of  only  80°C  when  the  sand/steel 
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interface  was  at  350°C.  It  then  continued  to  rise  very  slowly  for  some  time 
after  the  power  supply  was  removed.  The  very  low  temperature  on  the  upper 
surface  of  the  SCB  was  unexpected  and  so  we  ran  several  more  tests  to  try  to 
identify  the  cause  of  this  behaviour.  Figure  3b  is  a  repeat  run  with  the  SCB 
empty,  and  similar  results  were  obtained.  We  then  removed  the  screw-on  cap 
and  turned  the  SCB  upside  down.  Figure  3c  again  shows  the  temperature  on 
the  inside  of  the  steel  cylinder  (the  upper  curve),  while  the  initially  lower  curve 
is  the  temperature  recorded  at  the  centre  of  the  top  surface  of  the  upturned  SCB 
(i.e.  what  is  normally  the  base).  The  two  temperatures  are  now  much  closer 
and  *'..e  temperature  difference  between  the  two  points  is  due  to  the  finite  time 
t  ihen  for  the  heat  to  diffuse  from  the  heating  bands  to  the  upper  surface  (we 
describe  this  in  more  detail  when  discussing  Figure  7).  These  results  are 
important  because  they  show  that  the  large  mass  of  steel  in  the  cap  has  a 
significant  effect  on  the  heat  flow  and  temperature  distribution  within  the  SCB, 
and  indicate  that  the  effect  of  the  massive  steel  witness  plates  should  also  be 
included  in  any  numerical  simulation  of  the  full  SCB  assembly. 


Time,  seconds 


Figure  3a:  Temperature  measured  at  two  locations  within  an  SCB  filled  with  sand. 
Upper  curve  -  temperature  at  the  sand/steel  interface  at  a  position  halfway  up  the 
height  of  the  cylinder.  Lower  curve  —  temperature  on  the  top  surface  of  the  SCB  cap 
at  the  approximate  centre  point. 


13 


Time,  seconds 


Figure  3b:  As  for  Figure  3a,  except  that  the  sand  has  been  removed  and  the  SCB  is 
simply  air-filled. 


Figure  3c:  Temperatures  measured  at  two  locations  on  an  empty,  inverted  SCB  with 
the  screw  cap  removed.  The  (initially)  upper  curve  denotes  the  temperature  on  the 
inside  surface  of  the  steel  cylinder,  while  the  (initially)  lower  curve  is  the  temperature 
recorded  at  the  centre  of  the  top  surface  of  the  upturned  SCB. 


We  have  also  found  similar  behaviour  for  the  solid  steel  cylindrical  block 
described  in  the  previous  section.  The  SCB  heating  bands  were  placed  around 
the  cylinder  (in  exactly  the  same  position  as  when  applied  to  the  SCB,  refer  to 
Figure  1)  and  connected  to  a  240  volt  power  supply  for  120  seconds.  The  solid 
upper  and  lower  curves  in  Figure  4  are  the  temperatures  recorded  on  the 
curved  surface  of  the  block  midway  between  the  two  heating  bands  and  the 
temperature  at  the  centre  of  the  top  surface  of  the  block,  respectively.  The 
dashed  upper  and  lower  curves  are  the  temperatures  recorded  at  the  same 
locations  under  identical  experimental  conditions  for  another  run.  The 
difference  between  the  solid  and  dashed  curves  illustrates  the  degree  of 
reproducibility  between  the  two  runs.  The  SCB  heating  bands  do  not  cover  the 
entire  curved  surface  of  the  cylinder;  there  is  an  8  mm  gap  between  the  top 
edge  of  the  band  and  the  top  surface  of  the  block.  The  broken  upper  and 
lower  curves  in  Figure  4  are  the  temperatures  recorded  again  at  the  same 
locations,  only  this  time  tire  top  surface  of  the  upper  heating  band  has  been 
moved  up  to  be  flush  with  the  top  surface  of  the  block.  In  this  case  the  side 
temperature  is  lower  and  the  surface  temperature  higher  (after  a  certain  time) 
because  of  the  relative  placement  of  the  heating  band.  These  results  indicate 
that  very  accurate  numerical  simulations  of  the  temperature  distribution  within 
the  SCB  will  require  further  modification  to  the  HEAT  code  to  allow  the  finite 
size  of  the  heating  bands  to  be  specified. 


Figure  4:  Temperatures  measured  at  several  locations  on  a  solid  steel  cylinder  as  a 
function  of  time.  The  cylinder  has  the  same  dimensions  as  the  SCB  and  is  heated 
using  the  SCB  heating  bands.  The  solid  upper  and  lower  curves  are  the  temperatures 
recorded  on  the  curved  surface  of  the  block  midway  between  the  two  heating  bands,  and 
the  temperature  at  the  centre  of  the  top  surface  of  the  block,  respectively.  The  dashed 
upper  and  lower  curves  are  the  temperatures  recorded  at  the  same  locations  for  a 
replicate  run.  The  broken  upper  and  lower  curves  are  the  temperatures  recorded  again 
at  the  same  locations,  only  this  time  the  top  surface  of  the  upper  heating  band  has  been 
moved  up  to  be  flush  with  the  top  surface  of  the  block. 
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4.  HEAT  Code  Calculations 


Before  using  HEAT  to  model  the  complexity  of  the  multi-material  SCB  we  first 
applied  the  code  to  the  simpler  problem  of  modelling  the  temperature 
distribution  within  the  solid  cylindrical  steel  block.  HEAT  requires  as  input  a 
prescribed  time-dependent  temperature  profile  along  the  outer  curved  surface. 
For  the  solid  steel  block  we  used  an  average  of  the  two  uppermost  curves 
shown  in  Figure  4.  These  data  are  denoted  by  the  upper  solid  curve  in 
Figure  5.  The  filled  circles  denote  experimental  values  for  the  temperature  on 
the  surface  and  are  averages  of  the  two  lowermost  curves  in  Figure  4,  while  the 
lower  solid  curve  is  the  HEAT  code  result.  The  computational  grid  consisted  of 
125  cells  in  the  axial  direction  and  35  cells  in  the  radial  direction,  with 
Az  =  Ar  =  1  mm,  and  Af  =  10  s.  The  density  and  specific  heat  were  held 
constant  at  7860  kg/m3  and  490  J/kg.K  respectively,  and  the  thermal 
conductivity  was  varied  to  give  the  best  possible  fit  to  the  experimental  points. 
This  resulted  in  a  value  for  X  of  40.0  W/m.K.  Thermal  conductivity  varies 
greatly  between  steels  of  different  types,  and  the  value  determined  here  is 
consistent  with  values  quoted  in  the  literature  [7,  81. 


Figure  5:  Comparison  of  calculated,  and  measured  temperatures  for  the  solid 
cylindrical  block.  The  upper  solid  curve  is  the  average  measured  temperature  on  the 
side  of  the  block  and  is  used  as  the  input  surface  temperature  for  the  HEAT  code 
simulation.  The  filled  circles  denote  temperatures  measured  at  the  centre  of  the  top 
surface  of  the  block,  and  the  lower  solid  curve  is  the  HEAT  code  simulation.  The 
dashed  line  represents  temperatures  at  the  centre  of  the  top  surface  calculated  from  a 
one-dimensional  code. 
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Although  there  appears  to  be  reasonable  agreement  between  the  predicted 
and  experimental  temperature-time  curves  shown  in  Figure  5,  closer 
examination  shows  that  the  two  curves  have  quite  different  shapes,  indicating 
that  some  of  the  assumptions  used  in  the  numerical  model  may  not  be  valid. 
One  such  assumption  is  that  the  thermal  properties  of  the  materials  are 
temperature  independent.  This  is  discussed  further  in  the  next  section.  The 
approximate  agreement  between  the  two  curves  does  allow  us  to  conclude 
however  that  the  temperature  difference  between  the  top  and  side  surfaces  of 
the  SCB  observed  in  Figure  5  is  simply  due  to  the  finite  time  needed  for  the 
heat  to  diffuse  from  the  heating  bands  to  the  upper  surface  of  the  block. 

Experimentation  with  the  code  showed  that  the  results  were  insensitive  to  the 
values  used  for  e  and  h,  implying  that  heat  losses  through  the  ends  of  the 
cylinder  were  not  responsible  for  determining  the  temperature  distribution  over 
the  top  and  bottom  surfaces.  This  was  checked  by  calculating  the  surface 
temperature  using  a  simple  one-dimensional  explicit  code,  and  the  result  is 
shown  by  the  dashed  line  in  Figure  5.  The  temperatures  calculated  using  the 
one-dimensional  code  are  slightly  higher  than  those  calculated  using  HEAT,  but 
the  difference  is  at  most  a  few  degrees  C°.  This  clearly  shows  that  the  heat 
flow  along  the  axis  of  the  cylinder  has  little  effect  on  the  calculated 
temperatures,  and  demonstrates  that  the  primary  influence  on  the  surface 
temperature  profile  is  the  energy  input  from  the  heating  bands  on  the  surface  of 
the  cylinder,  and  their  location  from  the  point  at  which  the  temperature  is 
measured.  The  appropriate  HEAT  code  boundary  conditions  over  the  top  and 
bottom  surfaces  for  this  simulation  are  evidently 


3  T/dn  =  0, 


(5) 


rather  than  equation  (2),  although  the  changes  to  the  coding  to  implement 
equation  (5)  are  as  involved  as  those  required  to  implement  equation  (2).  It 
should  also  be  noted  that  this  particular  simulation  was  conducted  at  a  very 
rapid  heating  rate,  and  while  the  radiative  and  convective  heat  loss  terms  are 
insignificant  for  this  particular  heating  regime,  this  will  not  necessarily  be  the 
case  for  the  very  much  slower  heating  rates  used  for  slow  cookoff  studies, 
where  the  explosive  events  can  take  on  the  order  of  20  hours  to  occur. 

We  now  consider  the  temperature  within  an  inert  filled  SCB.  Figure  6  shows 
the  temperature  measured  at  the  sand/steel  interface  (the  initially  upper  solid 
line),  the  temperature  measured  at  the  centre  of  the  sand  fill  (the  filled  circles), 
and  the  HEAT  code  calculation  for  the  temperature  at  the  centre  of  the  sand 
using  two  different  values  for  the  thermal  conductivity  of  the  sand  (lower  solid 
line,  X  =  0.25  W/m.K;  dashed  line  X  =  0.27  W/m.K).  The  values  of  the 
physical  constants  used  for  the  three  materials  (sand,  steel  and  air)  are  listed  in 
Table  1.  The  discontinuity  in  the  slope  of  the  curve  for  the  temperature 
measured  in  the  centre  of  the  sand  is  caused  by  the  vapourization  of  adsorbed 
water  from  the  sand  at  100°C.  This  phase  change  effect  is  not  modelled  by 
HEAT  and  the  computed  curve  cannot  reproduce  this  discontinuity.  Provided 
comparisons  between  computed  and  experimental  curves  are  made  above  this 
temperature  however  this  will  not  affect  the  results  to  be  discussed  here. 
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Figure  6:  Comparison  of  calculated  and  measured  temperatures  for  the  inert  filled 
SCB.  The  (initially)  solid  upper  curve  is  the  temperature  measured  at  the  sand/steel 
interface  which  is  used  as  the  input  temperature  for  HEAT.  The  filled  circles  denote 
the  measured  temperature  at  the  centre  of  the  sand,  and  the  (initially)  lower  solid  curve 
and  the  dashed  curve  represent  the  HEAT  code  results  calculated  using  \  (sand)  equal 
to  0.25  and  0.27  W/m.K  respectively. 


The  agreement  between  calculated  and  experimental  temperatures  in  Figure  6 
is  reasonable,  and  if  the  errors  in  the  experimental  points  are  considered  to 
have  the  same  magnitude  as  those  shown  in  Figure  4  then  the  agreement  is 
very  close,  but  not  quite  exact.  The  two  curves  calculated  using  slightly 
different  values  of  thermal  conductivity  show  that  it  is  possible  to  improve 
agreement  in  some  parts  of  the  temperature-time  domain,  but  this  then  leads  to 
worse  agreement  in  other  parts.  It  was  found  impossible  to  vary  any  of  the 
material  parameters  in  such  a  way  as  to  improve  the  overall  agreement.  We 
attribute  the  inability  of  HEAT  to  predict  the  temperature  across  the  entire  time 
period  to  the  neglect  of  the  temperature  variation  of  the  thermal  parameters. 

We  subsequently  measured  the  specific  heat  of  the  sand  in  this  particular  test 
using  the  differential  scanning  calorimetry  method  and  found  it  to  vary  by  20% 
over  the  temperature  range  100°C  to  200°C  (see  Table  2).  We  anticipate  that  if 
HEAT  were  modified  to  include  the  temperature  dependence  of  the  thermal 
conductivity  and  heat  capacity  then  the  agreement  shown  in  Figure  6  would 
be  considerably  improved. 
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Table  1:  Material  Parameters  for  Inert  Filled  SCB 


Table  2: 


Sand: 

X 

=  0.25  W/m.K 

cP 

=  1108  J/kg.K 

P 

=  1410  kg/m3 

Ain 

X 

=  0.03  W/m.K 

=  1004  J/kg.K 

p 

=  1.2  kg/m3 

Steel: 

X 

=  40.0  W/m.K 

Cf 

=  490  J/kg.K 

p 

=  7860  kg/m3 

Heat  of  Inert  (Sand) 

Temperature  (K) 

C„  (J/kg.K) 

370 

825 

380 

843 

420 

896 

430 

917 

460 

944 

470 

972 

We  have  also  calculated  the  temperature  at  the  centre  of  the  sand  using  the 
one-dimensional  explicit  code  and  found  exact  agreement  with  the 
two-dimensional  HEAT  results.  This  illustrates  that  the  heat  loss  through  the 
end  surfaces  of  the  cylinder  has  little  effect  on  the  temperature  at  this  location 
and  for  this  heating  rate,  and  that  the  primary  influence  on  the  temperature  is 
provided  by  the  energy  flowing  into  the  sand  from  the  heating  bands. 

One  final  illustration  of  the  above  two  points  is  provided  by  Figure  7.  This 
shows  a  comparison  between  the  calculated  and  predicted  temperatures  for  the 
SCB  for  the  situation  described  in  relation  to  Figure  3c,  i.e.  the  SCB  is  empty 
and  inverted.  The  initially  upper  solid  curve  is  the  temperature  on  the  inside 
wall,  the  filled  circles  the  temperature  measured  on  the  top  of  the  base  at  the 
centre,  and  the  initially  lower  solid  curve  the  HEAT  code  prediction.  Again 
the  overall  agreement  between  calculated  and  measured  temperatures  is  only 
fair;  there  is  exact  agreement  for  early  times,  and  again  for  late  times  (although 
this  may  be  fortuitous),  but  it  was  found  impossible  to  vary  the  thermal 
parameters  within  realistic  limits  and  produce  agreement  everywhere. 
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Figure  7:  Comparison  of  calculated  and  measured  temperatures  for  the  SCB  for  the 
configuration  described  in  Figure  3c.  The  (initially)  upper  solid  curve  is  the 
temperature  on  the  inside  wall  of  the  SCB,  the  filled  circles  are  the  measured 
temperatures  on  the  top  of  the  SCB,  and  the  (initially)  lower  solid  curve  is  the  HEAT 
code  simulation  of  this  temperature.  The  dashed  curve  is  the  temperature  at  the  same 
point  calculated  using  a  one-dimensional  code. 

The  slightly  larger  discrepancy  between  the  calculated  and  experimental 
points  in  the  midpoint  of  the  range  may  also  be  due  to  the  effect  of  convective 
heat  flow  in  the  air  contained  within  the  SCB,  This  undoubtedly  occurs  in  the 
early  stages  of  heating  but  has  not  been  included  in  the  simulation,  which  may 
have  led  to  an  overall  shift  in  the  position  of  the  computed  curve  and 
exaggerated  the  discrepancy  at  the  midpoint.  The  temperature  has  also  been 
calculated  using  the  one-dimensional  code  and  this  is  shown  by  the  dashed  line. 
Again,  the  results  are  very  similar. 


5.  Discussion  and  Conclusion 

The  experimental  results  and  numerical  simulations  described  in  the  previous 
section  allow  us  to  make  the  following  observations: 


(a)  Accurate  numerical  simulation  of  the  temperature  distribution  within  some 
parts  of  the  SCB  will  not  be  possible  with  the  existing  version  of  HEAT. 

We  are  referring  in  particular  here  to  the  temperature  distribution  on  the 
surface  of  the  screw-on  cap  shown  in  Figures  3a  and  3b.  A  more  accurate 
prediction  could  be  made  by  modifying  parts  of  HEAT  to  allow  the  driving 
temperature  to  be  applied  only  to  those  parts  of  the  finite  difference  grid 
which  are  in  contact  with  the  heating  bands,  and  then  allowing  the 
remainder  of  the  curved  surface  to  lose  heat  through  radiation  and 
convection.  This  would  require  extensive  changes  however,  and  may  not 
be  necessary  for  our  purposes. 

(b)  The  temperature  within  the  explosive  sample  in  an  SCB  may  well  be 
determined  to  within  an  acceptable  accuracy  by  neglecting  heat  flow  along 
the  axis  of  the  SCB  and  using  a  one-dimensional  code  to  calculate 
temperature  along  the  radial  direction.  This  could  be  investigated  further 
using  FEAT  to  calculate  axial  temperature  profiles,  and  also  with 
additional  thermocouple  measurements  at  selected  points  within  the 
explosive  substitute. 

(c)  Accurate  prediction  of  the  temperature-time  profiles  within  the  explosive 
will  require  a  computer  code  which  allows  for  the  detailed  temperature 
dependence  of  the  explosive's  thermal  properties. 

Bearing  the  above  observations  in  mind,  we  make  the  following 

recommendations: 

(i)  Several  more  runs  should  be  made  on  the  inert  filled  SCB  using  additional 
thermocouples  to  obtain  a  better  knowledge  of  the  temperature  distribution 
within  the  explosive. 

(ii)  Additional  HEAT  code  simulations  should  be  performed  to  investigate  the 
temperature  variation  in  the  axial  direction. 

(iii)  The  temperature  dependence  of  the  thermal  parameters  for  the  inert  fill 
should  be  determined  accurately. 

(iv)  A  one-dimensional  heat  flow  code  which  takes  into  account  the 
temperature  dependence  of  the  material  properties  should  be  written  and 
then  applied  with  the  parameters  determined  in  (iii)  above  to  see  if 
improved  agreement  between  calculated  and  experimental  temperatures 
can  be  obtained. 

(v)  The  origin  of  the  instabilities  which  occur  for  some  combinations  of  Ar,  A z 
and  A/  should  be  investigated  more  thoroughly.  This  could  most  easily  be 
performed  by  considering  the  stability  of  the  finite  difference  scheme  across 
an  interface  using  a  one-dimensional  version  of  the  code. 
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